/*  sdfagmpi.cpp  by K.Tsuru  */
// function ID 3552     DRADIX only
#ifndef SN_H
#include "sn.h"
#endif
/*******************************************************
It Provides pi by Arithmetic and Geometric mean method
********************************************************/
static SDouble Sqrt2(const SDouble& a, const SDouble& b){
    RealSize C;
    uint sz = a.MaxSize();
    SDouble w(a.Type(), sz), y(a.Type(), sz), one(1.0);
    static SDouble recsqrt(0.0);
    uint ef = (DOUBLE_FIG*2u)/DFIGURES, fig = a.EffFig();
    int c = DDCompare(a, b);
    if(!c) return a;    // a == b

    long agfig = AgreeUpTo(a, b);
    w = a*b;
    if(agfig < (long)ef*DFIGURES) return Sqrt(w);
    ef = (uint)max((long)ef, (2*agfig)/DFIGURES);
    ef = min(ef, fig);

    int power = w.RdxExp()*DFIGURES/2;
    w.MultPowRdx( -w.RdxExp() );

/* set an initial value y0 */
    if( w < 0.01 ) { w *= 100.0;  power--; }
    C.SetEffFig(ef);
    if(recsqrt.Sign()){
        y = recsqrt;
    }else {
        double w0=doubleD(w);
        y = 1./sqrt(w0);
    }
    C.SetEffFig(0);

    SDouble delta(a.Type(), sz);
    w.FixedPoint(w.RdxExp());

    int fp = 0, itrmax = howpow2( (DFIGURES*y.MaxSize())/DOUBLE_FIG+1 ) +6;
    c = 0;
    do{
        if( (ef = C.SetEffFig(ef) ) >= fig) fp = 1;
        delta = one-w*(y*y);
        delta *= y;
        delta /= 2;
        y += delta;
        c++;
        ef *= 2;
        if(ef > fig) ef = fig;
        C.SetEffFig(0);
    } while( (!delta.IsMLT(y) && (c < itrmax)) || !fp );

    recsqrt = y;
    y.MultPow10(-power);
    y = w*y;
#ifndef NDEBUG
    if(c == itrmax) y.SetError(y.FATAL, "Sqrt2()", 3552);
#endif
//  y.iterationCount = c;
    y.PointFree();
    y.Reform(10011);

    return y;
}
/// Using Gauss-Legendre method ///
SDouble AGMPi(){
    SDouble a(1.0), b(2.0), t(0.25), x(1.0), y, d, pi;
    int c;
    b = RecSqrt(b);
    y.FixedPoint(0);
    int itrmax = 2*howpow2(a.MaxSize());
    c = 0;
    do{
        y = a;
        a = (a+b)/2.0;
        b = Sqrt2(b, y);
//      fprintf(stderr, "Iteration counts = %2d\n", b.ItrCounts());
        d = y-a;
        d *= d;
        t = t - x*d;
        x *= 2;
        d = a - b;
        c++;
    } while(!d.IsMLT(a) && (c < itrmax) );
#ifndef NDEBUG
    if(c == itrmax) y.SetError(y.FATAL, "AGMPi()", 3552);
#endif
    d = a+b;
    d *= d;     t *= 4;
    pi = d/t;
    pi.iterationCount = c;
    return pi;
}